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We consider the problem of state selection for a stochastic system, initially in an unstable 
stationary state, when multiple metastable states compete for occupation. Using path-integral 
techniques we derive remarkably simple and accurate formulas for state-selection probabilities. The 
method is sufhciently general that it is applicable to a wide variety of problems. 
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The investigation of the decay from a metastable state 
has been the subject of numerous studies over very many 
years [0-^ . But the analogous problem of the decay from 
an unstable state has received comparatively little atten- 
tion |^|-^, and what studies there have been have fo- 
cussed on the kinetic properties of one-dimensional and 
quasi-one-dimensional systems [Q. However, these stud- 
ies cannot address one of the fundamental, open ques- 
tions in non-equilibrium statistical mechanics: state se- 
lection from an unstable state in systems with multiple, 
isolated minima. Here we present a systematic, intuitive, 
and analytically tractable method which gives results in 
excellent agreement with Monte-Carlo simulations. 

When driven far from equilibrium many systems en- 
counter instabilities. At such points, noise plays a crucial 
role. In addition, in complex systems there are multiple 
modes that interact and can compete. Perhaps the most 
familiar example is found in Rayleigh-Benard convection. 
Consideration of the interaction between two competing 
modes leads to the following equations for their ampli- 
tudes X and y B] : 
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where a and f3 are the (positive) growth rates for the two 
modes x and y, 7 is the (positive) coupling coefficient, 
and S and e are positive stabilizing coefficients. The vari- 
ables rjx and rjy are Gaussian random variables with mean 
zero and variance {'ni{t)rjj{t')) = 2D6ijS{t — t'), where i 
and j are either x 01 y, and D is the noise strength. As 
we are considering the decay from the unstable stationary 
point x — Q,y = Q, the noise plays an essential role. There 
are four main elements that are present in Eq. (Q): (i) an 
unstable stationary point with exponential growth of the 
modes in the neighborhood of this point, (ii) interaction 
between the modes, (iii) isolated metastable states, and 
(iv) noise. These features are also found in many other 
systems 

In this Letter we address the question: given a sys- 
tem described by equations such as (Q), with the initial 
condition being the unstable stationary point, what is 
the probability that the system finds itself in a given 
metastable configuration? The system will relax to ther- 



mal equilibrium over a time scale that is on the order of 
exp{E / D) , where is a characteristic energy barrier sep- 
arating the metastable states. However, if ^ 1, this 
time can be enormous. Our focus is on understanding 
the occupation over shorter time scales. 

Our approach is based on the path-integral represen- 
tation for the conditional probability density P{r, T|0, 0) 
that the system resides in state r at time T given that 
it started at the unstable stationary point at the ori- 
gin. For purposes of illustration, we take the concrete, 
physically important example presented in Eq. ([T]). The 
path-integral expression for P is given by [021 
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Here S is the action, J is the Jacobian, and V is the 
potential for this problem and is given by 
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To evaluate the path integral for weak noise, a natural 
approximation scheme is the method of steepest descent. 
In this approach, the path integral is dominated by the 
paths of least action; a necessary condition is that these 
paths make the action stationary. The leading approxi- 
mation is to simply evaluate the action along these paths. 
However, in our case it is necessary to go beyond this or- 
der and include both (Gaussian) ffuctuations about the 
paths of least action, as well as the Jacobian evaluated 
along the appropriate path. In other words, once the 
stationary paths have been determined, we need to cal- 
culate three quantities: (i) the action, (ii) the Jacobian, 
and (iii) the fluctuation determinant, which characterizes 
the effect of ffuctuations about the relevant path. 
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it end at = {X, 0) at i = T. The analogous problem 
where the endpoint is Vy = (0, Y) can be handled in ex- 
actly the same way.) The most obvious stationary path 
is Tc = {xc,0), where 
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The action Sc — S[xc\ for this solution is given by 
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FIG. 1. (a) V{x,y) with a = l3 = 2, ■y = 4, S = e= 1/5. 
(b) Contours of zero force. 



In the limit that T — > cxd, 5 — > 0, so that, at least in 
this limit, this solution is a path of least action, not 
simply a stationary path. We are unable to prove that 
this is the only stationary path that connects (0, 0) with 
{X,0). However, we will show that we can make sig- 
nificant progress by considering only this path. In fact, 
this simplification will enable us to derive a remarkably 
simple formula for P{rx, T\0, 0). 

The second factor that we must evaluate is the Jaco- 
bian evaluated along the path Vc- Using Eqs. (jj) and (Q) 
we find that 



Having outlined a general prescription for calculating 
the conditional probability density, we now focus on the 
specific example introduced above. Figure la is a three- 
dimensional plot of V(r) for a certain choice of param- 
eters. In Fig. lb we plot the locus of points for which 
VV^ = 0. This is a useful way to visualize state space; 
the points where the ellipses intersect each other are the 
saddle-points of V (denoted by crosses) , the points where 
the ellipses intersect the x and y axes are the local minima 
(denoted by open circles), and the origin is the unstable 
stationary point (denoted by a closed circle). The ques- 
tion of interest here can now be phrased in the following 
way: given an ensemble of systems, each of which starts 
at the unstable stationary point, what fraction of the 
ensemble flows into an x-valley or y- valley (which lead 
to the X- and y- wells, respectively)? For simplicity, we 
suppose that 6 and e are sufficiently small that the local 
minima are so distant from the region where state selec- 
tion occurs that they have no influence. Operationally, 
this consists of setting S and e to zero, so that now 
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Given that the minima are now irrelevant, it is natural 
to consider the conditional probability that if the sys- 
tem starts at r = it ends up in an x-valley denoted 
by {X, 0) or a y- valley denoted by (0, Y). To do this we 
follow the procedure outlined above, viz. we first find the 
path, or paths, of least action that connect the unstable 
point to a point in one of these valleys. (Hereafter we 
shall confine our attention to the calculating the proba- 
bility that if the system starts at r = (0, 0) at t = that 
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It is straightforward to calculate both and Jj,, 
the results that 



Jx = exp 
and Jy = exp 
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The third quantity that we must calculate is the ef- 
fect of fluctuations about Vc- To do this, we expand 5'[r] 
about the path Tc, keeping terms of second order. Taking 
V = Vc + i5r, we have that ^[r] = S[vc\ + ^ / dtSr L[yc] (5r, 
where Lc = L[rc] is a 2 x 2 matrix-differential operator 
that is given by 
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The path integral in Eq. (|2|) over r now becomes an in- 
tegral over (Sr. Using the second-order expansion of S'[r], 
the Gaussian integrals over 5y can be completed. These 
integrals contribute a factor of y^jdetLd to the ex- 
pression for the conditional probability. As L is block- 
diagonal, we have that detLc ~ dci LxAei Ly. Com- 
bining this with the fact that Jc = JxJy, the steepest- 
descent approximation for the conditional probability can 
be written as 
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The expressions for Sc, Jx and Jy are given in Eqs. {jg), 
( |lO| ) and (^l|), respectively. The conditional probability, 
P, is the product of two factors: Pq, which is described 
in the following paragraph, and fl^^. This form is par- 
ticularly appealing because, as we shall see below, it is f2 
that accounts for the presence of the competing y-mode, 
whereas Pq is independent of both (3 and 7 and is there- 
fore insensitive to the presence of y. To determine Pq 
and fl we need to calculate det and det Ly . 

The calculation of det is straightforward, with the 
result that detL^, oc sinh(aT). Combining this result 
with Eq. (^) for Sc and Eq. ( p^ ) for J^, Pq can be writ- 
ten as 
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Pq is the conditional probability density that a one- 
dimensional system under the influence of the potential 
—ax'^/2 and Gaussian white noise be located a,t x = X 
at t = T given that it started at cc = at i = 0. As a 
function of T, Pq is peaked at a value T* that is given 
by coth(ar*) = 1 + 2D/{aX'^), i.e. 



T* = {2a)-^ ln(l + aX^/D). 
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The calculation of det Ly is not as straightforward; it 
can be expressed as [n3| 
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where hi and /i2 are two linearly independent solutions 
of the homogeneous equation Lyh = 0. The denomina- 
tor of Eq. (gO|) is the Wronskian of the two solutions. To 
evaluate Eq. (^), consider the quantity 

hi{X,t) = exp (^pt-j dt' Xc{t'f^ . (21) 

Taking the second derivative of hi with respect to t we 
find that 
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Comparing Eq. ( |2^ ) with Lyh = from Eq. ([i4|), we 
see that if Xc = axc then hi is one of the desired so- 
lutions. Now Xc = aXcCoth{at), so as long as coth(ai) 
is close to 1, /ii is a good solution. Recall, however, 
that for T < T*, Po is essentially zero [c.f. Eq. (^]. 
Thus, it is only values of T > T* that are relevant. But 
coth(aT*) = 1 -|- 0{D) so indeed we expect that as long 



as D <C 1, hi is a good approximate solution to the 
homogenous equation Lyh = 0. The second linearly in- 
dependent solution, /i2 can be expressed in terms of hi 
as h2{t) = hi{t) Jp dt' hi'^{t'). With this choice of /12, we 
have that /i2(0) = 0. In addition, hi{0) = 1 so that the 
Wronskian is unity and 
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Jo 



(23) 



Combining this equation with the fact that Jy — 
h1^'^{T) [c.f. Eqs. (P) and (|2^)], we find that 
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With this expression for fi, together with Eq. (|T|) for 
hi{t) and Eq. ( p^ for Pq, we have succeeded in deriving 
a formula for P{rx, T\0, 0) that accounts for the Jacobian 
prefactor as well as the Gaussian fluctuations about the 
stationary path Vc- To calculate the analogous formula 
for P{ry, T\0, 0) [where Vy = (0, Y)], we simply switch a 
and f3 and replace X with Y. 

At this stage we are in a position to calculate, us- 
ing our expression for P{rx,T\O,0), the probability that 
the system flows into a given well. Currently, we have 
a simple analytic expression for Pq, but to calculate 
O we need to integrate h^^ over time. This inte- 
gration presents no difficulty in principle. However, 
by making two approximations we are able to obtain 
a simple analytic formula for Q. Specifically, we ap- 
proximate the exponential factor exp[^^^^ — sinh(2at) 
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1 + exp 

tribution from the lower fimit (T = 0) in Eq. (|2^). This 
second approximation is due to the fact that, as explained 
in the discussion following Eq. (|2^), hi(t) is only a good 
solution for at > 1. We now obtain the following approx- 
imate formula for il{X,T): 
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We have compared Eq. (|25| ) with numerical calculations 
of Eq. (HJ), and we find that for the relevant range of 
parameters the results are essentially indistinguishable. 
Thus, Eq. ([l|) for Pq and Eq. (H) for n together pro- 
vide an analytic formula for the conditional probability 
given in Eq. (p^. 

We now turn to the computation of the relative prob- 
ability that the system flows into an x- or y-valley. The 
strategy is to calculate the total probability flux through 
the cc-valleys and the y-valleys and compare them. The 
probability current, which we denote by J^{r, T), is given 
by J = -PVy-P)VP, so that the total flux through 
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an X- valley at X is Tx{X) = J^dt 
the total flux J^y through a y-valley at 



^dyX(r,t) and 
Y is Ty{Y) = 
/o I-ao^y'^vi^^^)j where J'x and J'y are the x- and 
y-components of J, respectively. Denoting by the 
relative probability of flowing into an x-valley we then 
have that 
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The calculation of the fluxes requires a knowledge of 
J and hence of P(r, r|0,0). In particular, we require 
this function for an arbitrary point in the x-valley and 
not just on the x-axis, i.e., we require -P(r), not simply 
P{vx)- The method we have presented may be extended 
to obtain the full functional dependence on r but 
the results given above do not give P any explicit de- 
pendence on the y variable across the x-valley or the x 
variable across the j/- valley. We will therefore limit our- 
selves here to showing that, by estimating the flux by 
sampling it on the axis, we can get excellent agreement 
with Monte-Carlo simulations and therefore conflrm the 
essential correctness of our approach. A feature of this 
procedure is the necessity of fitting X or Y . We expect 
that this will no longer be required when the flux is cal- 
culated, since this should enable J-x{X) and J-y{Y) to 
be calculated for large X and Y where we would expect 
them to be insensitive to their actual values. The only 
restriction that we will impose on X and Y is that they 
are not too small, for then state-selection will not have 
occurred when these points are reached. We estimate the 
minimum value of X to be of the order of Xmin, the point 
at which the force in the y-direction changes sign. For 

V given in Eq. (||) we have that Xmin = \/WH- Like- 
wise, we have that the minimum value of y is given by 
Ymin = \foLpi. 

The Monte-Carlo simulations are performed on the 
Langevin equation with V{x, y) given in Eq. (^. In Fig. |^ 
the results are shown for a range of values of a and par- 
ticular choices of f3 and of (7 and D always appear in 
this combination, since the effect of the interaction is to 
renormalize the noise) . The theory we have outlined here 
is seen to be in excellent agreement with the simulations. 
For = 0.1 we have taken X = Xmin and Y = Ymin 
and for 7D = 0.001 we have taken X = 1.83Xmin and 

Y = l.83Ymin- Comparison for other values of the pa- 
rameters, a determination of the region of validity of our 
approximation in parameter space and further improve- 
ments of the method will also be discussed elsewhere . 

In this Letter we have presented a systematic method 
for determining state selection from an unstable station- 
ary state, when multiple metastable states compete for 
occupation. Previous methods have not addressed this 
question directly. Our treatment has the added advan- 
tage of yielding closed form, analytic expressions for the 



conditional probability distribution. Finally, we empha- 
size that, although we have focussed on a specific po- 
tential system with two degrees of freedom for illustra- 
tive purposes, our theory is neither restricted to potential 
problems nor to systems with only two degrees of free- 
dom. 
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FIG. 2. Probability of flowing into an x-valley as a func- 
tion of Q. Simulation results are for l3 — 1 and the continuous 
curves are our theoretical results. 

We thank Ken Elder for useful discussions, and the 
Universities of Chicago and Manchester for hospital- 
ity. This work was supported in part by EPSRC grant 
GR/K79307 (AJM) and by the NSF (DMR-9415604) 
(MBT). 



[8 

[9 
[10: 

[11 
[12: 

[is: 

[14 



H. A. Kramers, Physica (Utrecht) 7:284 (1940). 

H. Risken, The Fokker-Planck Equation (Springer- 

Verlag, Berlin, 1989), Section 5.10. 

C. W. Gardiner, Handbook of Stochastic Methods 

(Springer- Verlag, Berlin, 1985), Chapter 9. 

M. Suzuki, J. Stat. Phys. 16, 11 (1977). 

B. CaroU, C. CaroU and B. Roulet, J. Stat. Phys. 26, 83 

(1981). 

U. Weiss, Phys. Rev. A 25, 2444 (1982). 
By quasi-one-dimensional we mean the potential is spher- 
ically symmetric and the system is therefore effectively 
one- dimensional . 

L. Segel, J. Fluid Mech. 14, 97 (1962); A. NeweU and J. 
Whitehead, J. Fluid Mech. 38, 279 (1969); R. Graham, 
Phys. Rev. A 10, 1762 (1974). 

M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65 
851 (1993). 

L. Kramer, H. R. Schober and W. Zimmermann, Physica 
D 31, 212 (1988). 

B. Bagchi and G. R. Fleming, J. Phys. Chem. 94, 9 
(1990). 

R. Graham in Noise in Dynamical Systems Vol. 1, eds. F. 
Moss and P. V. E. McClintock (CUP, Cambridge, 1989). 
A. J. McKane and M. B. Tarlie, J. Phys. A 28 6931 
(1995). 

A. J. McKane and M. B. Tarlie. In preparation. 



4 



